### Set Directory
try(setwd(dirname(rstudioapi::getActiveDocumentContext()$path)))
rm(list=ls())

### Load Package
library(stargazer)
library(dplyr)
library(tidyr)

### Load Data
data = read.csv('replication dataset.csv')

# Restrict sample to potential targets
# that have economic sanctions experienced
# threat and/or imposition during 1975-2005
data_pf = data %>%
  filter(year %in% 1975:2005) %>%
  arrange(ccode,year) %>%
  group_by(ccode) %>%
  mutate(ever = ifelse(mean(threat)>0|mean(imposition)>0,1,0)) %>%
  ungroup() %>%
  filter(ever==1)

### First Stage OLS
m1 = lm(v2juncind.rev~v2juncind.rev.lag+hrdummy+nonhrdummy+polconiii+yr_sch_sec+open.log, data=data_pf)
m2 = lm(v2x_corr~v2x_corr.lag+hrdummy+nonhrdummy+polconiii+yr_sch_sec+open.log, data=data_pf)
m3 = lm(v2mecenefm.rev~v2mecenefm.rev.lag+hrdummy+nonhrdummy+polconiii+yr_sch_sec+open.log, data=data_pf)
# Table 3
stargazer(m1,m2,m3,type='text')

### Second Stage Bootstrap
set.seed(2020)
N=1000

# Function to get fitted value
hat = function(d,m){
  X = d %>%
    select(names(coef(m))[-1])
  X = cbind(1,X)
  fit = as.matrix(X)%*%coef(m)
  return(fit)
}

## PTS A

# leave out observations with any missing data
da = data_pf %>%
  select(v2juncind.rev.lag,v2x_corr.lag,v2mecenefm.rev.lag,
         hrdummy,nonhrdummy,polconiii,yr_sch_sec,open.log,
         tptsa,PTS_A,
         gdppclog,polity2,civilwar,interwar) %>%
  drop_na()

# Function to bootstrap second stage
boot1 = function(x){
  d1 = da[sample(1:nrow(da),replace=T),]
  df1 =d1 %>%
    mutate(fit1 = hat(d1,m1)) %>%
    mutate(fit2 = hat(d1,m2)) %>%
    mutate(fit3 = hat(d1,m3))
  
  mb1 = lm(tptsa~PTS_A+fit1+fit2+fit3
           +gdppclog+polity2+civilwar+interwar,
           data=df1)
  
  beta1 = coef(mb1)
  
  return(beta1)
}

# Bootstrapping (N=1000 iterations)
betas1 = replicate(N,boot1())

# 95% Confidence Interval
result1 = t(apply(betas1,1,function(x){
  coef = mean(x)
  lwr = sort(x)[floor(N*.025)]
  upr = sort(x)[floor(N*.975)]
  lwr2 = sort(x)[floor(N*.05)]
  upr2 = sort(x)[floor(N*.95)]
  return(cbind(coef,lwr,upr,lwr2,upr2))
}))

result1 = as.data.frame(result1)

colnames(result1) = c('coef','lwr','upr','lwr2','upr2')
result1$id = 1:nrow(result1)

# Table OA3 Model 8
round(result1[,1:3],3)
nrow(da)

# Figure 3 Left Panel
par(mfrow=c(1,2),
    mar=c(5,2,3,2))

plot(id~coef,data=result1[3:5,],ylim=c(6,2),
     xlim=c(0,.5),yaxt='n',xlab='Coefficient',ylab='',pch=c(15,16,17),
     main='PTS A')
abline(v=0,lty='dashed')
for(i in 3:5){
  lines(y=c(i,i),x=c(result1[i,2],result1[i,3]))
}
text(y=result1[3:5,6]-0.2,x=c(0,.2,0),labels=c('Lower Court Dependence','Political Corruption','Government Censorship'),pos=4)

## PTS S

# leave out observations with any missing data
ds = data_pf %>%
  select(v2juncind.rev.lag,v2x_corr.lag,v2mecenefm.rev.lag,
         hrdummy,nonhrdummy,polconiii,yr_sch_sec,open.log,
         tptss,PTS_S,
         gdppclog,polity2,civilwar,interwar) %>%
  drop_na()

# Function to bootstrap second stage
boot2 = function(x){
  d1 = ds[sample(1:nrow(ds),replace=T),]
  df1 =d1 %>%
    mutate(fit1 = hat(d1,m1)) %>%
    mutate(fit2 = hat(d1,m2)) %>%
    mutate(fit3 = hat(d1,m3))
  
  mb1 = lm(tptss~PTS_S+fit1+fit2+fit3
           +gdppclog+polity2+civilwar+interwar,
           data=df1)
  
  beta1 = coef(mb1)
  
  return(beta1)
}

# Bootstrapping (N=1000 iterations)
betas2 = replicate(N,boot2())

# 95% Confidence Interval
result2 = t(apply(betas2,1,function(x){
  coef = mean(x)
  lwr = sort(x)[floor(N*.025)]
  upr = sort(x)[floor(N*.975)]
  lwr2 = sort(x)[floor(N*.05)]
  upr2 = sort(x)[floor(N*.95)]
  return(cbind(coef,lwr,upr,lwr2,upr2))
}))

result2 = as.data.frame(result2)

colnames(result2) = c('coef','lwr','upr','lwr2','upr2')
result2$id = 1:nrow(result2)

# Table OA3 Model 9
round(result2[,1:3],3)
nrow(ds)

# Figure 3 Right Panel
plot(id~coef,data=result2[3:5,],ylim=c(6,2)
     ,xlim=c(0,.5),yaxt='n',xlab='Coefficient',ylab='',pch=c(15,16,17),
     main='PTS S')
abline(v=0,lty='dashed')
for(i in 3:5){
  lines(y=c(i,i),x=c(result2[i,2],result2[i,3]))
}
text(y=result2[3:5,6]-0.2,x=c(0,.2,0),labels=c('Lower Court Dependence','Political Corruption','Government Censorship'),pos=4)

## FIW PR

# leave out observations with any missing data
dpr = data_pf %>%
  select(v2juncind.rev.lag,v2x_corr.lag,v2mecenefm.rev.lag,
         hrdummy,nonhrdummy,polconiii,yr_sch_sec,open.log,
         tpr,pr,
         gdppclog,polity2,civilwar,interwar) %>%
  drop_na()

# Function to bootstrap second stage
boot3 = function(x){
  d1 = dpr[sample(1:nrow(dpr),replace=T),]
  df1 =d1 %>%
    mutate(fit1 = hat(d1,m1)) %>%
    mutate(fit2 = hat(d1,m2)) %>%
    mutate(fit3 = hat(d1,m3))
  
  mb1 = lm(tpr~pr+fit1+fit2+fit3
           +gdppclog+polity2+civilwar+interwar,
           data=df1)
  
  beta1 = coef(mb1)
  
  return(beta1)
}

# Bootstrapping (N=1000 iterations)
betas3 = replicate(N,boot3())

# 95% Confidence Interval
result3 = t(apply(betas3,1,function(x){
  coef = mean(x)
  lwr = sort(x)[floor(N*.025)]
  upr = sort(x)[floor(N*.975)]
  lwr2 = sort(x)[floor(N*.05)]
  upr2 = sort(x)[floor(N*.95)]
  return(cbind(coef,lwr,upr,lwr2,upr2))
}))

result3 = as.data.frame(result3)

colnames(result3) = c('coef','lwr','upr','lwr2','upr2')
result3$id = 1:nrow(result3)

# Table OA4 Model 10
round(result3[,1:3],3)
nrow(dpr)

# Figure OA1 Left Panel
plot(id~coef,data=result3[3:5,],ylim=c(6,2)
     ,xlim=c(0,.35),yaxt='n',xlab='Coefficient',ylab='',pch=c(15,16,17),
     main='FIW PR')
abline(v=0,lty='dashed')
for(i in 3:5){
  lines(y=c(i,i),x=c(result3[i,2],result3[i,3]))
}
text(y=result3[3:5,6]-0.2,x=c(0,.1,0),labels=c('Lower Court Dependence','Political Corruption','Government Censorship'),pos=4)

## FIW CL

# leave out observations with any missing data
dcl = data_pf %>%
  select(v2juncind.rev.lag,v2x_corr.lag,v2mecenefm.rev.lag,
         hrdummy,nonhrdummy,polconiii,yr_sch_sec,open.log,
         tcl,cl,
         gdppclog,polity2,civilwar,interwar) %>%
  drop_na()

# Function to bootstrap second stage
boot4 = function(x){
  d1 = dcl[sample(1:nrow(dcl),replace=T),]
  df1 =d1 %>%
    mutate(fit1 = hat(d1,m1)) %>%
    mutate(fit2 = hat(d1,m2)) %>%
    mutate(fit3 = hat(d1,m3))
  
  mb1 = lm(tcl~cl+fit1+fit2+fit3
           +gdppclog+polity2+civilwar+interwar,
           data=df1)
  
  beta1 = coef(mb1)
  
  return(beta1)
}

# Bootstrapping (N=1000 iterations)
betas4 = replicate(N,boot4())

# 95% Confidence Interval
result4 = t(apply(betas4,1,function(x){
  coef = mean(x)
  lwr = sort(x)[floor(N*.025)]
  upr = sort(x)[floor(N*.975)]
  lwr2 = sort(x)[floor(N*.05)]
  upr2 = sort(x)[floor(N*.95)]
  return(cbind(coef,lwr,upr,lwr2,upr2))
}))

result4 = as.data.frame(result4)

colnames(result4) = c('coef','lwr','upr','lwr2','upr2')
result4$id = 1:nrow(result4)

# Table OA4 Model 11
round(result4[,1:3],3)
nrow(dcl)

# Figure OA1 Right Panel
plot(id~coef,data=result4[3:5,],ylim=c(6,2)
     ,xlim=c(0,.35),yaxt='n',xlab='Coefficient',ylab='',pch=c(15,16,17),
     main='FIW CL')
abline(v=0,lty='dashed')
for(i in 3:5){
  lines(y=c(i,i),x=c(result4[i,2],result4[i,3]))
}
text(y=result4[3:5,6]-0.2,x=c(0,.1,0),labels=c('Lower Court Dependence','Political Corruption','Government Censorship'),pos=4)

### Summary Statistics

# Table OA2
stargazer(as.data.frame(data_pf%>%
                          select(PTS_A,PTS_S,pr,cl,
                                 v2juncind.rev,v2x_corr,v2mecenefm.rev,
                                 hrdummy,nonhrdummy,polconiii,yr_sch_sec,open.log,
                                 gdppclog,polity2,civilwar,interwar))
          ,type='text')
